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1.  Introduction 


Extensions  of  the  original  constant-temperature  Dissipative  Particle  Dynamics  (DPD)  method 
( 1 ,  2)  have  been  developed,  including  methods  that  impose  constant-pressure  (DPD-P)  (3,  4) 
conditions.  Formulated  from  an  extended  system,  the  DPD-P  method  can  be  implemented  using 
either  a  Hoover  or  Langevin  barostat  (2).  DPD-P  allows  control  over  both  the  temperature  and 
the  pressure,  where  in  typical  fashion  the  cell  volume  fluctuates  to  satisfy  the  imposed  pressure. 

When  applying  any  DPD  method,  numerical  integration  of  the  equations-of-motion  (EOM)  is  a 
key  consideration  since  the  stochastic  component  requires  special  attention.  Efficient  and 
accurate  integration  schemes  are  required  to  maintain  reasonable  time  step  sizes,  thus  allowing 
for  the  simulation  of  actual  mesoscale  events.  Moreover,  the  integration  is  a  nontrivial  task  due 
to  the  pairwise  coupling  of  particles  through  the  random  and  dissipative  forces  (5).  In  recent 
work  ( 6 ),  a  comprehensive  description  of  numerical  integration  schemes  based  upon  the 
Shardlow-splitting  algorithm  (SSA)  was  presented  for  the  constant-temperature  DPD  method. 
The  original  SSA  formulated  for  systems  of  equal-mass  particles  was  extended  to  systems  of 
unequal-mass  particles.  Both  a  velocity- Verlet  scheme  and  an  implicit  scheme  were  formulated 
for  the  integration  of  the  fluctuation-dissipation  contribution  where  the  velocity- Verlet  scheme 
consistently  performed  better. 

In  this  work,  we  formulate  the  SSA  for  the  DPD-P  method,  where  we  verify  the  method  using 
both  the  standard  DPD  fluid  (pure  and  binary  mixtures)  (7)  and  a  coarse-grain  solid  model  (5). 
For  completeness,  derivations  of  the  Fokker- Planck  equation  (FPE)  and  the  fluctuation- 
dissipation  theorem  (FDT)  are  included. 


2.  Formulations  of  DPD  at  Fixed  Temperature  and  Pressure  Using 
Shardlow-Like  Splitting  Numerical  Discretization 


2.1  General  Formulation  of  DPD 

DPD  particles  are  defined  by  a  mass  mj ,  position  r, ,  and  momentum  p, .  The  particles  interact 
with  each  other  via  a  pairwise  force  F„  that  is  written  as  the  sum  of  a  conservative  force  F~  , 
dissipative  force  F? ,  and  random  force  F*  : 

Fi=F?+k’+V  a) 

F(  is  given  as  the  negative  derivative  of  a  coarse-grain  potential,  //  6 ,  expressed  as 
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where  ry.  =  r(.  -  r;.  is  the  separation  vector  between  particle  i  and  particle  j  and  r;  =  |r;  | .  The 
remaining  two  forces,  F/|)  and  F*  ,  can  be  interpreted  as  a  means  to  compensate  for  the  degrees- 

of-freedom  neglected  by  coarse-graining.  The  conservative  force  is  not  specified  by  the  DPD 
formulation  and  can  be  chosen  to  include  any  forces  that  are  appropriate  for  a  given  application, 
including  multibody  interactions  (e.g.,  \8-10\).  Fy/;  and  FyK  are  defined  as 

/  A 


K  = 


-y a ® 


r 

v 


(3) 


and 


Fy  =  W 


r 

v 


(4) 


where  yu  and  aij  are  the  friction  coefficient  and  noise  amplitude  between  particle  i  and  particle 

j  ,  respectively,  v ..  =  —  -  —  ,  and  Wn  are  independent  Wiener  processes  such  that  Wu  =  WH  . 
m  m 

‘  J 

The  weight  functions  coD{r)  and  coR{r)  vanish  for  r>rc  where  rc  is  the  cut-off  radius. 

Note  that  F(>c  is  completely  independent  of  F;|)  and  F*  ,  while  F^  and  F^  are  not  independent 

but  rather  coupled  through  a  fluctuation-dissipation  relation.  This  coupling  arises  from  the 
requirement  that  in  the  thermodynamic  limit,  the  system  samples  the  corresponding  probability 
distribution. 


2.2  Constant-Pressure  DPD 


Generally,  a  barostat  in  an  extended  system  approach  is  introduced  via  variables  £ ,  p!: ,  and  WE . 
£  is  defined  as  the  logarithm  of  the  system  volume  V  ,  s  —  In  where  l7 (O)  is  the  volume  at 
t  -  0,  WE  is  a  mass  parameter  associated  with  £  ,  and  p£  is  the  momentum  conjugate  to  £  (11). 


A  Langevin  barostat  can  be  incorporated  into  this  extended  system  approach  via  additional 
dissipative  and  random  terms  (12).  The  DPD-P  method,  first  introduced  by  Jakobsen  (3),  was 
formulated  for  both  a  Hoover  and  Langevin  barostat,  where  soon  after  Trofimov  et  al.  presented 
a  similar  formulation  based  upon  an  Andersen  barostat  (4).  For  uniform  dilation  using  a 
Langevin  barostat,  the  EOM  are  given  as 
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instantaneous  pressure,  P0  is  the  imposed  pressure,  d  is  the  dimensionality  of  the  system, 


Nf  =  dN  —  d,YP  and  ap  are  the  Langevin  barostat  parameters,  and  Wp  is  the  Wiener  process 
associated  with  the  random  fluctuations  of  the  piston.  Note  that  when  yp  -  ap  =  0,  the  EOM 
shown  in  equation  5  reduce  to  the  EOM  corresponding  to  the  extended  system  approach  (i.e.,  a 
Hoover  barostat)  (13).  W£  is  usually  expressed  as  W£  =  [n f  +d)k^T t2p  ,  where  zp  is  the 


characteristic  time  of  the  barostat  that  should  be  chosen  slightly  larger  than  the  smallest  time 
scale  of  the  particle  motions  (14).  Analogous  to  constant-temperature  DPD,  the  dissipative  and 
random  force  parameters  must  conform  to  the  FDT  corresponding  to  the  isothermal-isobaric 
probability  density.  These  constraints  are  satisfied  by  imposing  the  constant-temperature  FDT 
relations 


<j2  =  2yijkBT 

^  (6) 

coD{r)=\coR{r)l 

along  with  a  FDT  relating  the  piston  parameters,  given  as 

a2P=2ypW£kRT,  (7) 


where  kB  is  the  Boltzmann  constant  and  T  is  the  temperature.  The  FPE  and  an  outline  of  the 
derivation  of  the  FDT  are  presented  in  appendix  A.  With  further  analogy  to  constant-temperature 
DPD,  the  conservation  of  the  total  momentum  is  again  due  to  pairwise  additivity  of  the 
conservative,  dissipative,  and  random  forces.  Furthermore,  in  the  limit  of  yt]  — »  0  and  yp  — »  0, 
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the  system  conserves  the  quantity  E'—  K  +  U  +  P0V  -+ - —  where  K  =  ^  is  the  kinetic 

2W£  V  2/77,. 

energy  and  U  =  IZ  is  the  configurational  energy.  For  completeness,  the  EOM  for 

i  j>i 

nonuniform  dilation  is  given  in  appendix  B. 
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2.3  Numerical  Discretization 

Jakobsen  ( 3 )  integrated  the  EOM  given  in  equation  5  using  a  modified  velocity- Verlet  algorithm. 
In  this  work,  we  propose  the  numerical  discretization  of  equation  5  using  a  splitting  strategy 
similar  to  that  employed  for  the  constant-temperature  DPD  formulation  presented  previously,  the 
Shardlow- splitting  velocity- Verlet  algorithm  (SSA-VV)  (6).  The  deterministic  differential 
equations  and  the  elementary  stochastic  differential  equations  (SDEs)  corresponding  to  equation 
5  are  as  follows.  The  conservative  terms  are 
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while  the  fluctuation-dissipation  terms  have  the  same  expressions  as  the  constant-temperature 
DPD  formulation, 
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where  the  superscript  i  -  j  indicates  that  the  variation  of  momenta  is  considered  for  a  pair  of 
interacting  particles  i  and  j  only,  and  d Wtj  -  d Wn  are  the  increments  of  the  Wiener  processes. 


The  stochastic  flow  map  <f>At  can  be  approximated  by  (6,  15) 
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where  °  denotes  a  composition  of  operators.  For  each  term,  momenta  updates  are  based 
upon  the  constant-temperature  DPD  formulation  previously  given  (<5): 
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where  the  superscript  i  -  j  has  been  omitted  for  notational  simplicity,  gtj  -  c/;  is  a  Gaussian 


random  number  with  zero  mean  and  unit  variance  that  is  chosen  independently  for  each  pair  of 
interacting  particles,  and  jufj  =  —  +  — . 


m,.  m; 


<f>At  can  be  treated  using  the  velocity- Verlet  scheme  proposed  by  Martyna  et  al.  (13, 14),  which 

requires  the  calculation  of  quantities  at  both  t  +  At  and  t  +  — ,  followed  by  an  iterative  process 

to  determine  the  remaining  quantities  at  t  +  At .  This  scheme  proceeds  by  first  solving  the 
following  set  of  equations: 
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where  L(o)  and  L(l  +  At)  are  the  lengths  of  the  cubic  simulation  box  at  t  -  0  and  t  +  At, 
respectively.  Next,  the  conservative  forces  at  t  +  At,  {f/  (7  +  At))^=l ,  are  evaluated  and 
subsequently  used  in  the  second  part  of  the  algorithm,  which  requires  an  iterative  approach.  The 
iteration  starts  with  an  estimation  of  pe  at  t  +  At  using  pf\t  +  At)  =  p£  (t  -  At)  +  2 A tFs (?) , 
followed  by  solving  the  set  of  equations 


){k\t  +  At)=- 


exp  [s(t  +  At)  -  a{t)  p; 


(  A  ^ 
t  +  — 

V  2  j 


1  + 


At 


f  -  'n  (*-i 


2  + 


d 


V  NfJ 


+  ^{t  +  At) 

-  (i  =  l,...,N) 


Pe  \t  +  kt) 

VK 


Fjk\t  +  At)  =  dv(t  +  At)[p^k\t  +  At)- P0]+ 
-rPPe~l)(t  +  At)+ffp-J= 

Pe\t  +  =Ps(t  +  +  +  At) 

V  2)2 


d  vpf)(t  + At)-p (k\t  +  At) 


(12b) 


m. 


self-consistently  until 


_ i _ 

3N 


is  less  than  a  prescribed  tolerance,  which 


is  typically  less  than  10  6 .  In  equation  (12b),  gp  is  a  Gaussian  random  number  with  zero  mean 
and  unit  variance  and  ( k )  is  the  iteration  index. 


The  practical  implementation  of  the  SSA-VV  approach  for  the  DPD-P  variant  is  analogous  to  the 
constant-temperature  DPD  formulation  ( 6 )  with  the  exception  that  the  deterministic  integration 
steps  are  replaced  by  equations  12a  and  12b.  For  completeness,  the  numerical  discretization  for 
nonuniform  dilation  is  given  in  appendix  B. 


1.  Stochastic  Integration  for  all  i  -  j  pairs  of  particles 
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1  +  Mij  y.ofAt  L 

^  •  i.l 


TJL.y. 

r  * 
Vv  J 

—  •v 


r... 


V  v  ) 


r.  2  r. 

ij  v 


—+^rcrij0}R£ij  —  ^ 

ra  2  n, 


+  cr,jcoRgl  T,J 


a  yfAt 


ru  2 


r  r  Vat 

-cr.-OJ  cu— - 

r,  2 


2.  Deterministic  Integration  #1.1 


(i)  P: 


f 


At 


t  + 

V  2y 


=  Pe(t)+ 
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(ii)  <?(r  +  Ar)  =  +  Ar  - 


dpt 


At 

t  H - 

V _ 2  j 

W 


(iii)  V(t  +  At)  -  V(o)exp  \s(t  +  A/1)] 

(iv)  L(t  +  At  )  —  V(t  +  At  j 


V/3 


3.  Deterministic  Integration  #1.2  for  i  =  1,...,  N 


f 


(i)  P, 


At 


t  H - 

v  2  j 


t  \  At 

=  P,W  +  y 


Pit)- 


2  + 


d_ 

Nr 


w. 


PPM 


(ii)  r  (t  +  At)  =  exp \s(t  +  At)- s(t) 


){t)+At- 


P, 


At 
t  +  — 

V _ 2y 

m; 


4.  Conservative  Force  Calculation:  {f/ 

5.  Deterministic  Integration  #2.1 

pf\t  +  At)  =  p,(t-  A  t)+  2  AtF.:{t) 


6.  Deterministic  Integration  #2.2  for  i  = 

a)  pi\t+At)= 


exp  [e(t  +  At) -  £-(r)]pj  it  +  ^  F,c  (t  +  At) 

V  2)2 


1+ 


At 


2  + 


d_ 

Nt 


'pFC  At) 


FE(k\t  +  At)  =  clV(t  +  Ar)[pu)(r  +  At)- P0]+ 


d  vpf)(/1  +  At)-p(k)(t  +  At) 


(ii) 


/  >' 


m, 


■YpP{s~'](t  +  At)  +  (jp 


Va ~t 
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7.  Deterministic  Integration  #2.3 

At 


r 


t  H - 

v  2  J 


At 


-F^\t  +  At) 


Ps\t  +  At)=p 
Steps  6  and  7  are  carried  out  self-consistently  until 

prescribed  tolerance. 


Z[p!‘y+A,)-p!‘-"(,+A,)f 


3  N 


is  less  than  a 


3.  Computational  Details 


The  SSA-VV  for  the  DPD-P  method  was  tested  using  both  the  standard  DPD  fluid  (7)  and  a 
coarse-grain  solid  model  (<8);  complete  details  of  the  conservative  forces  for  these  models  are 
given  in  appendix  C.  Both  a  pure  component  case  and  an  equimolar  binary  mixture  were  tested 
for  the  DPD  fluid  model.  System  sizes  for  the  DPD  fluids  and  coarse-grain  solid  were, 
respectively,  N  =  10125  and  13500.  For  these  simulations,  the  following  reduced  units  were 
used:  rc  and  r0  are  the  unit  of  length  for  the  DPD  fluid  and  coarse-grain  solid,  respectively;  the 
mass  of  a  DPD  particle  is  the  unit  of  mass;  and  the  unit  of  energy  is  kBT  .  Using  these  reduced 
units,  we  set  the  maximum  repulsion  between  particles  i  and  j  as  atj  =  25  for  the  pure  DPD 
fluid.  For  the  binary  DPD  fluid,  the  values  were  a{]  -  25  and  28  for  the  like  and  unlike  i  -  j 
interactions,  respectively.  Further,  for  all  cases,  we  set  the  noise  amplitude  cr  =  3  and  the 
barostat  characteristic  time  tp  -  2.  Prescriptions  for  the  choice  of  yp  (3,  14)  suggest  that  the 
value  should  be  between  2!  rp  and  10/  rp ;  therefore,  we  set  yp  =  10/ rp  =  5  for  all  cases. 


4.  Results 


As  a  test  of  the  SSA-VV  formulation,  we  verify  that  the  DPD-P  variant  converges  to  the  same 
equilibrium  properties  when  at  the  same  thermodynamic  conditions  as  a  constant-temperature 
DPD  simulation. 

4.1  DPD  Fluid 

The  benchmark  systems  for  both  the  pure  and  binary  DPD  fluid  cases  are  taken  from  a  constant- 
temperature  DPD  simulation  performed  at  p  =  3  and  T  =  1  and  run  for  trun  =  3000  and 
At  =  0.03 .  The  following  quantities  were  evaluated  and  are  listed  in  table  1  (pure  fluid)  and 
table  2  (equimolar  binary  fluid):  virial  pressure  (Pvir) ,  configurational  energy  per  particle  (uj , 
kinetic  temperature  (Tkin) ,  and  self-diffusion  coefficients  I)  using  the  Einstein  relation  (11). 
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To  validate  the  constant-pressure  SSA-VV  formulations,  DPD-P  simulations  were  performed 

using  both  the  Langevin  and  Hoover  barostats  with  an  imposed  pressure  determined  from  the 
constant-temperature  DPD  simulation  (P0  =  23.65  and  PQ  =  24.79  for  the  pure  and  equimolar 

binary  DPD  fluids,  respectively)  for  trun  =  3000  and  A t  =  0.03 .  The  results  for  the  pure  and 
equimolar  binary  DPD  fluids  are  summarized  in  tables  1  and  2,  respectively,  where  excellent 
agreement  between  the  DPD-P  and  constant-temperature  DPD  simulations  are  found  for  (u) , 
(Tkin } ,  the  particle  density  (p"'j  =  >  ar|d  D  .  Asa  further  test,  DPD-P  simulations  were  also 

started  from  a  random  configuration  subject  to  an  energy  minimization  (as  opposed  to  starting 

from  the  equilibrated  configuration  of  the  constant-temperature  DPD  simulation),  where  again 
calculated  quantities  were  in  excellent  agreement  with  constant-temperature  DPD  results  (not 
shown). 

Table  1.  The  configurational  energy  per  particle  (llj  ,  the  kinetic  temperature  (Tkinj , 
the  virial  pressure  ( Pvjr  'j  ,  the  particle  density  (^p'j  ,  and  the  self-diffusion 
coefficient  D  ,  determined  from  simulations  of  the  pure  DPD  fluid.  (}j 
denotes  an  ensemble  average,  where  numbers  in  parentheses  are  uncertainties 
calculated  from  block  averages. 


Variant 

(«) 

(TJ 

(p*) 

(P) 

D 

DPD 

P  =  3 

4.56(1) 

1.005(8) 

23.65(8) 

— 

0.295(13) 

DPD-P  Langevin 

P0  =  23.65 

4.55(1) 

1.005(8) 

23.59(8) 

2.997(8) 

0.294(14) 

DPD-P  Hoover 

P0  =  23.65 

4.55(1) 

1.004(8) 

23.61(8) 

2.997(8) 

0.296(19) 

Table  2.  The  configurational  energy  per  particle  (ll^ ,  the  kinetic  temperature  {^kin  )  ,  the 
virial  pressure  (Pv!r)  ,  the  particle  density  (^p'j  ,  and  the  self-diffusion 
coefficient  D  ,  determined  from  simulations  of  the  equimolar  binary  DPD  fluid. 
Q  denotes  an  ensemble  average,  where  numbers  in  parentheses  are 
uncertainties  calculated  from  block  averages. 


Variant 

(») 

(O 

<o 

(p) 

D, 

d2 

DPD 

P  =  3 

4.76(1) 

1.005(8) 

24.79(13) 

— 

0.177(13) 

0.165(13) 

DPD-P  Langevin 
P0  =  24.79 

4.76(2) 

1.004(8) 

24.76(4) 

2.998(8) 

0.176(15) 

0.163(17) 

DPD-P  Hoover 

P0  =  24.79 

4.76(1) 

1.004(8) 

24.75(4) 

2.998(9) 

0.179(12) 

0.165(9) 
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4.2  Coarse- Grain  Solid 


A  validation  study  analogous  to  the  DPD  fluids  study  above  is  performed  for  a  coarse-grain  solid 
model,  where  a  recently  developed  nickel  model  is  considered  that  reasonably  reproduces  several 
measured  properties,  including  the  melting  temperature  (8).  For  a  benchmark  system,  a  constant- 
temperature  DPD  simulation  is  performed  at  p  =  8260  kg/m  and  T  =  1300  K  for  trun  =  1  ns  and 
At  =  5  fs,  where  results  are  listed  in  table  3.  (Since  the  coarse-grain  solid  model  has  been 
parameterized  to  an  actual  material,  results  are  reported  in  real  units  as  opposed  to  reduced  units 
for  the  DPD  fluid.)  At  this  state  point,  the  atomistic  Sutton-Chen  model  of  nickel  predicts  a 
pressure  of  approximately  0  bar  (8),  while  (Pvir)  for  the  coarse-grain  solid  model  is  larger  than  0 
bar.  Starting  from  an  equilibrated  configuration  from  a  constant-temperature  DPD  simulation, 
nonuniform  dilation  DPD-P  simulations  were  performed  using  both  the  Langevin  and  Hoover 
barostats  at  P0  =  0  bar  for  trun  =  1  ns  and  At  -5  fs,  where  results  are  given  in  table  3.  Compared 
to  the  constant-temperature  DPD  case,  the  DPD-P  results  are  in  near  exact  agreement  for  both 
barostats. 


Table  3.  The  molar  configurational  energy  (ll'j ,  the  kinetic  temperature  \Tknt )  , 
the  virial  pressure  ( Pvjr  \ ,  and  the  mass  density  / p \  ,  determined  from 
simulations  of  the  coarse-grain  solid  model  of  nickel,  (.'j  denotes  an 
ensemble  average,  where  numbers  in  parentheses  are  uncertainties 
calculated  from  block  averages. 


Variant 

(«) 

(kj/mol) 

(T») 

(K) 

(P+) 

(bar) 

(p) 

(kg/m3) 

DPD 

p  =  8260  kg/m3 

-508.44(11) 

1300.1(91) 

5.91(94) 

— 

DPD-P  Langevin 

P0=  0  bar 

-508.43(14) 

1299.9(92) 

-0.13(95) 

8259.3(73) 

DPD-P  Hoover 

P0  =  0  bar 

-508.44(14) 

1299.7(91) 

0.06(99) 

8260.1(68) 

5.  Conclusion 


Comprehensive  descriptions  of  numerical  integration  schemes  based  upon  the  SSA  were 
presented  for  the  isothermal-isobaric  DPD  method.  The  SSA  was  readily  extendable  to  both  the 
Hoover  and  Langevin  barostats  under  both  uniform  and  nonuniform  dilation,  where  the  SSA  for 
all  variants  was  found  to  be  a  stable  and  accurate  integration  scheme.  The  equivalence  of  the 
DPD  variants  was  verified  using  both  a  standard  DPD  fluid  model  and  a  coarse-grain  solid 
model,  where  thermodynamic  quantities  were  considered. 
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Appendix  A.  Fokker-Planck  Equation  (FPE)  and  Fluctuation-Dissipation 

Theorem  (FDT) 


The  FPE  corresponding  to  the  equations  of  motion  (EOM)  given  by  equation  5  of  the  report  is 


Lcp  +  LDp  +  LLBp  , 
ot 


(A-l) 


where  the  conservative  operator  is 


-Jmi  dr,.  ,J  dp,.  J  W£  i  dr, 


~  —  T>,-  •  —  -  1  +  —  —  YPt 

X\7  t—i  1  O-  AT  W 


d  P 


N  f  )W£  i  dp, 


dV  _d_ 
d t  8V 


mi  dp 


•  (A-2) 


The  operator  that  accounts  for  the  effects  of  the  dissipative  and  random  forces,  LD ,  is  given  by 


i  j*i  1  ij  ^r/ 


r> ico  + 


a  r, 


dp,,  dp,.  rt 


(A- 3) 


where  mi  is  the  mass  of  particle  i ,  r  =  r,  —  r  .  is  the  separation  vector  between  particle  i  and 
particle  j  ru  =  |ry  |  ^  is  the  conservative  force  acting  between  particle  i  and  particle  J  , 

5  9 

Y ij  and  <Tij  are  the  friction  coefficient  and  noise  amplitude  between  particle  i  and  particle  j 

respectively,  v,  =  —  -  —  ,  and  (0D  and  coR  are  weight  functions  of  the  dissipative  and  random 
,J  m,  /  7i  j 

forces,  respectively.  The  operator  representing  the  Langevin  barostat  terms  in  the  EOM  is 


^LB  — 


YPPs  + 


a],  d 
2  dp£ 


(A-4) 


In  equations  (A- l-A-4):  P  =  p(r»p»  V  is  the  system  volume,  Pe  is  the  momentum 

conjugate  to  £  =  In >  V(o)  is  the  volume  at  t  =  0,  W£  is  a  mass  parameter  associated  with 

s ,  P  is  the  instantaneous  pressure,  P0  is  the  imposed  pressure,  d  is  the  dimensionality  of  the 
system,  Nf  =  dN  -  d  _  and  Yp  and  a p  are  Langevin  barostat  parameters. 
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In  addition  to  the  system’s  implicit  contact  with  a  heat  reservoir  as  in  the  previous  constant- 
temperature  case,  constant -pressure  dissipative  particle  dynamics  introduces  a  barostat  that  keeps 
the  system  pressure  constant  through  implicit  contact  with  a  piston.  Under  these  circumstances, 
the  equilibrium  probability  density  peq  =  peq  (r,  p,  V,  p,: )  then  corresponds  to  the  isothermal- 
isobaric  probability  density 


P"1  (r .  p,  V ,  p£ )  =  ~  exp 


„2 

tf(r,p)+P0V  +  i- 

KT 


=  — exp 
Q 


Pl  ' 


ZU^+ZZ  “F+w  +  w 

i  2mt  ij?i  _ 2^ 

kj 


(A- 5) 


where  Q  is  the  normalizing  partition  function.  Similarly,  as  before,  Lcpeq  =  0  since  the 
isothermal-isobaric  probability  density  is  the  equilibrium  solution  for  the  conservative  system. 
The  FDT  then  follows  from  the  requirements  that  LDpeq  =  0  and  LLBpeq  =  0  .  The  former  leads 
to1 

rl 

r 

while  the  latter  is  satisfied  for 

a;  =  2ypWekJ .  (A-7) 


2  LjkJ 

coD 


(A-6) 


Note  that  the  temperature  of  the  Langevin  barostat  corresponds  to  the  temperature  of  the  heat 
reservoir,  which  maintains  the  system  temperature.  Therefore,  the  amplitude  of  the  volume 
fluctuations  depends  not  only  on  the  Langevin  barostat  parameters  but  also  on  the  system 
temperature. 


1  Brennan,  J.  K.;  Lfsal,  M.  Dissipative  Particle  Dynamics  at  Isothermal  Conditions  Using  Shardlow-Like  Splitting 
Algorithms ;  ARL-TR-6582;  U.S.  Army  Research  Laboratory:  Aberdeen  Proving  Ground,  September  2013. 


15 


Intentionally  lelt  blank. 
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Appendix  B.  Constant  Pressure  Dissipative  Particle  Dynamics  (DPD-P)  for 

Nonuniform  Dilation 


For  constant-temperature  and  constant-pressure  nonuniform  dilation  using  a  Langevin  barostat, 
the  equations  of  motion  are  given  as 


*ia  =  ^Jt  +  ^riadt 

l,a  nr  W,  l,a 


dp  =Y(fc  +Fd  +Fr  W/t-^p  dt 

Jl  i,a  /  j\  ij,a  ij,a  ij,a  J  -ywr  Jri,a 

j*i  We 

Ps,x+Ps.y+Ps 


£,x  e,y  a  e,z 

3NfWe 


Pi,adt 


d\nh=^dt 
WE 

^Ps.a  =  Kadt 


(i  =  1,...,N)  , 


(B-l) 


where:  a  =  (x,  y,z),  sa  =  In — j—r ,  La  and  La( 0)  are  the  box  lengths  for  the  a  -direction  at  t 
and  t  =  0,  respectively,  p..  a  is  the  momentum  conjugate  to  sa ,  Ws  is  the  mass  parameter 
associated  with  sa ,  typically  set  to  We  =  “(Y  +  v2p  ,  rp  is  the  barostat  characteristic  time, 

N,=3N-3, 


P  =I 

°P  y 


z 


PuaPuP 


V  »  mi 


■YYFc  r 

ij,a  ij.[ 


~ ~  L  - -  -  rPPe,a  +  °PWP,a  > 

Nf  i  mi 


is  the  component  of  the  pressure  tensor  for  the  a  and  /? 


i  j>‘  j 

directions,  P0  is  the  imposed  pressure,  V  -  LxLyL, ,  yp  and  ap  are  the  Langevin  barostat 
parameters,  aP  =  2ypW£kBT  is  the  FDT  relating  the  piston  parameters,  and  WP  a  is  the  Wiener 
process  associated  with  the  random  fluctuations  of  the  piston  in  the  a  -direction. 

B.l  Numerical  Discretization 

Applying  a  numerical  integration  splitting  strategy  similar  to  the  uniform  dilation  DPD-P  variant, 
the  deterministic  differential  equations  and  the  elementary  stochastic  differential  equations 
(SDEs)  corresponding  to  equation  B-l  are  as  follows.  The  conservative  terms  are 


Pi,a 


ir  =^il  +  Eyr  d, 

m ,  vv„ 


dt 


Ps.a  Pe,x  + Pe,y+ Pe,z 


J*l 


dlnL  =—dt 


3  NfWe 


Pi,adt  (i  =  l...,N) 


(B-2) 


dPe,a  =  Fs,adt 
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while  the  fluctuation-dissipation  terms  consist  of  the  elementary  SDEs: 


dp"7 

dp'-7 

j 


f 


-YiP 


l.y 

IJ 

r. 

V  v  J 


"  d/  -  CJ.  C)K  ”  dW 


RVU 


-dP' 


f  \ 

for  each  i  <  j 

V  J 


(B-3) 


As  previously  stated,  the  stochastic  flow  map  <f>At  is  approximated  by  equation  10  of  the  report. 
For  each  fluctuation-dissipation  term  ,  momenta  are  updated  by  applying  the  same 
expressions  as  in  the  uniform  dilation  case,  i.e.,  equations  1  la,  1  lb,  1  lc,  and  1  Id  of  the  report. 


tf>Al  can  be  treated  using  the  velocity- Verlet  scheme  proposed  by  Martyna  et  al.,>  which  requires 

the  calculation  of  quantities  at  both  t  +  At  and  t  +  ,  followed  by  an  iterative  process  to 

determine  the  remaining  quantities  at  t  +  At .  This  scheme  proceeds  by  first  solving  the  following 
set  of  equations 


A 

t  H - 

2  j 


=  Ps,Jf)+^FsJf) 


M  +  At)-sa(t)+  At- 


t  +  - 


At 


w„ 


La{t  + At)  =  La  (o)exp  [sa  (t  +  At)] 

V  {t  +  A  t )  —  Lx  (t  +  A  t  )Lv  ( [t  +  A  t  )Lz  (t  +  A  t ) 


,  At 
Pi, a  t  + 


=  Pi,a(t)+  y 


pCM  ^Ps,a(t)n  (A  P<Jt)+  A,v(0+  Pjt)  n  (A 

Fi.a  V )  ~  2  Pi, a  V) - tTAF - Pi,a  V) 


w„ 


3  NfWt 


ri,a  (*  +  At)  =  exp  [sa  (t  +  At)  -sa(t) 


Pi.a\  t  + 


At 


m. 


.  (B-4a) 


Next,  the  conservative  forces  at  t  +  At ,  |F(c  (/  +  A t)}*=l ,  are  evaluated  and  subsequently  used  in  the 
second  part  of  the  algorithm,  which  requires  iteration.  The  iteration  process  starts  with  an 
estimation  of  p£  a  at  t  + At  using  p  fl  (t  +  At)  —  p£  a  (t  -  At)  +  2 A tFe  a  (t) ,  followed  by  solving  the 

set  of  equations 


1  Martyna.  G.  J.;  Tobias,  D.  J.;  Klein,  M.  L.  Constant  Pressure  Molecular  Dynamics  Algorithms.  J.  Chem.  Phys.  1994,  101, 
4177. 

2 

~  Martyna,  G.  J.;  Tuckerman,  M.  E.;  Tobias,  D.  J.;  Klein,  M.  L.  Explicit  Reversible  Integrators  for  Extended  Systems 
Dynamics.  Mol  Phys.  1996,  87,  1117-1157. 
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/®*+A0=- 


exp  k  (t  +  a/ )  - -  (Ok«  (j  +  y  ] +  y  F£  (f  +  A?) 

,  |  AtL/kk+A0 ,  pl*;1V+A0+pl!;1)(f+A0+pl*;1)(f+A0' 


3NfWs 


F${t  +  At)  =  V{t  +  At\p^(t  +  A  0  -  d  +  TT  Z  P^  +  Af)'P^  +  Af) 


(B-4b) 


■^pk)(f+AO+^- 


(f  +  At )  =  P,,«  f ?  +  +  T  F<S (f  +  Af ) 

V  l  )  l 


self-consistently  until  — - — -  is  less  than  a  prescribed  tolerance,  which 

is  typically  chosen  to  be  0{  10~6 ).  In  equation  (B-4b),  gPa  is  a  Gaussian  random  number  with 
zero  mean  and  unit  variance  and  ( k )  is  the  iteration  index. 
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For  the  models  considered  in  this  work,  the  details  of  the  conservative  forces  expressed  in 
equation  2  of  the  main  text  are  the  following.  uGG  for  the  pure  and  binary  Dissipative  Particle 

Dynamics  (DPD)  fluids1  is  given  by 

UyG  =aijrcO)D(rij)  ,  (C-l) 

where  at/  is  the  maximum  repulsion  between  particle  i  and  particle  j  . 

For  the  coarse-grain  solid  model,  which  has  a  face-centered-cubic  (f.c.c.)  lattice  structure, 
particles  interact  through  a  shifted-force  Sutton-Chen  embedded  potential  (SC)  given  as 


where 


CG 

u;  =8 


j5ti 


f  X 

XX-  *  . 

y ) 

Pi  =  LPij  > 


P»  =  w(r»)-w(r,i  =  rt)-^7^U  (ru~r.)  • 


w 


(  \ 

(r)  = 

r 

~ 

V  «/ 

r. 

'v  y  J 

(C-2) 


^root,  R.  D.;  Warren,  P.  B.  Dissipative  Particle  Dynamics:  Bridging  the  Gap  Between  Atomistic  and  Mesoscopic 
Simulation.  J.  Chem.  Phys.  1997.  107,  4423. 
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where  s  and  r0  are  the  energy  and  length  parameters,  respectively,  n  and  m  are  positive 
integers  (n  >  m  to  satisfy  elastic  stability  of  the  crystal),  and  c  is  a  dimensionless  parameter. 
Although  effectively  this  is  a  many-body  potential,  the  force  on  each  particle  can  be  written  as  a 
sum  of  pairwise  contributions.  The  coarse-grain  solid  model  used  here  approximates  nickel  (Ni), 
where  one  DPD  particle  was  chosen  to  represent  four  f.c.c.  unit  cells,  i.e.,  16  Ni  atoms.  SC 
potential  parameters  were  determined  by  fitting  to  various  0-K  properties  and  the  melting 
temperature  at  zero  pressure,2  where  the  following  values  were  found:  si kB  -  225  K,  rQ  =  8.8698 
A,  c-  39.4314,  m  —  6 ,  and  n  =  9  .  Further  details  for  determining  SC  parameters  based  upon 
such  a  procedure  can  be  found  elsewhere.2'3 


2 

“Brennan,  J.  K.;  Lisal,  M.  Coarse-Grain  Models  for  Metals:  Constant-Pressure  Dissipative  Dynamics  Simulations.  In 
Proceedings  of  the  14th  International  Detonation  Symposium ,  Coeur  d’Alene,  ID,  1 1-16  April  2010;  Office  of  Naval  Research, 
2010;  p  1451. 

2Sutton,  A.  P.;  Chen,  J.  Long-Range  Finnis-Sinclair  Potentials.  J.  Philos.  Mag.  Lett.  1990,  61  (3),  139. 
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Intentionally  lelt  blank. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


DPD 

DPD-P 

EOM 

f.c.c. 

FDT 

FPE 

SC 

SDE 

SSA 

SSA-VV 


constant-temperature  Dissipative  Particle  Dynamics 

constant-temperature,  constant-pressure  Dissipative  Particle  Dynamics 

equations  of  motion 

face-centered-cubic 

fluctuation-dissipation  theorem 

Fokker-Planck  equation 

Sutton-Chen  embedded  potential 

stochastic  differential  equation 

Shardlow- splitting  algorithm 

Shardlow- splitting  algorithm-velocity  Verlet 


25 


NO.  OF 

COPIES  ORGANIZATION 


1  DEFENSE  TECHNICAL 
(PDF)  INFORMATION  CTR 
DTIC  OCA 

1  DIRECTOR 

(PDF)  US  ARMY  RESEARCH  LAB 
IMAL  HRA 

1  DIRECTOR 

(PDF)  US  ARMY  RESEARCH  LAB 
RDRL  CIO  LL 

1  GOVT  PRINTG  OFC 
(PDF)  AMALHOTRA 

ABERDEEN  PROVING  GROUND 


1  DIR  USARL 
(PDF)  RDRL  WML  B 
J  BRENNAN 


26 


